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A PROGRAM FOR CONTOURING RANDOMLY SPACED DATA 


By Roy W. Hamm, James F. Kibler, and W. Douglas Morris 
Langley Rosearch Center 

SUMMARY 

This paper presents a brief description of a digital computer program 
which prepares contour plots o'f three-dimensional data. As presently 

nfigured, the program can accept up to 56,000 randomly spaced data points, 
though the required computer resources may be prohibitive. However, with 
relatively minor internal modifications, the program can handle essentially 
unlimited amounts of data. The contouring techniques use a triangulation 
procedure developed by Dr. C, L. Lawson of the Jet Propulsion Laboratory. 

The contour points can be fitted with a smooth curve using an interpolating 
spline under tension. Sample cases illustrate contours of ocean wave 
simulation data. A general description of the main program and primary level 
subroutines is included to permit simple modifications of the program to 
handle other applications. 

INTRODUCTION 

Large amounts of three-dimensional data have been generated as a result 
of parametric studies of wave conditions off the Mid-Atlantic Coast (Ref. 1) . 
One way to analyze the data is by contouring wave height and wave orbital 
velocity parameters at selected values. In this way, the areas of differing 
wave activities can be defined. This information, in turn, is valuable in 
locating future offshore facilities in areas of minimum wave activity, or for 



locating dumping sites for either maximum or minimum dispersion of the dumped 
effluent or spoil. However, because of the large number of data points 
(80,000 generated by the study) and since these points do not fit into a 
i*ectangular grid pattern, existing Langley Research Center (LRC) contouring 
programs cannot adequately handle the data. 

In order to satisfy these requirements, a digital computer program has 
been modified to contour wave data. The program was originally developed by 
Dr. C. L. Lawson of the Jet Propulsion Laboratory. His "Triangulation" 
technique allowed the contouring of randomly spaced input data, without first 
■fitting the data into a rectangular grid (required my most other programs) . 
Later modifications were made by Dr. Alan K. Cline of ICASE*. One of the 
modifications incorporated a spline under tension to fit a smooth curve to 
the contour points. T'us program was then adapted to the Langley Research 
Center CDC 6600 system. Further modifications have been made which reduce 
the memory storage, execution time, and increase the plotting efficiency. 

The modified .contouring program represents an improvement over models 
presently available at LRC because it can provide smooth contouring of a 
relatively large number of randomly spaced data points (56,000). With minor 
internal modifications, the program could handle larger amounts of data. 

This paper presents a brief description of the contouring technique, a 
description of the computer program, and sample results from four typical 
cases . 


^Institute for Computer Applications in Science and Engineering, NASA-Langley 
Research Center 
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Description of Contouring Technique 


A contour plot is a convenient way to display 
a function of two independent variables. The con- 
tour is simply a line connecting constant values 
of the function plotted on a grid composed of the 
corresponding independent variables. In this way, 
a three-dimensional surface can be easily visual- 
ized in two dimensions. This section briefly 
describes the triangulation contouring technique. 

Assume that the function to be contoured is 

available as an array of randomly spaced data 

points. Each data point consists of an X,Y and 

Z value. The X and Y values correspond to the 

independent coordinates to be plotted and Z = 

Z(X,Y) is the function value to be contoured. The 
contouring process begins by forming a convex region 

composed of triangles from the set of X-Y data 
points Cfig* I)- Figure 2 illustrates an arrange- 
ment of triangles which form a convex region from 
the data points in figure 1. Obviously, this 
arrangement is not unique. In 7 rder to resolve 
this problem, a criterion is established to maxi- 
mize the smallest interior angle between possible 
pairs of triangles. In addition, the criterion is 


Y 


• • • 

X 

Figure 1.- A set of input data 
points . 


Y 



Figure 2.- A convex region 
composed of 
triangles . 



required to more evenly distribute interpolated contour points. For example, 
given a set of four data points, there are two choices for arranging the two 
triangles (figure 3). Under the above criterion, 3-b is chosen as the preferred 
arrangement . 




Figure 3.- Two choices for arranging triangles for four data points. 


After the triangles are arranged, each Z value is assigned to its 
corresponding X-Y vertex in each triangle. An interpolating plane which 
fits tho vortices exactly is then calculated for each triangle. If a 
chosen contour level lies between the Z values for 
any triangle, the two. contour points which lie on 
the sides of the triangle are computed from the 
interpolating plane (fig. 4) . These contour points 
are calculated for each triangle in the convex 
region and for each specified contour level. The i5,8 
contour points for any given contour level may then 
be connected by straight lines to yield a polygonal 



approximation to the required contour. Alternatively. ... 

Figure 4.- A triangle with 

a cubic spline function under tension (ref. 2) can linearly inter- 

polated contour 

i points. 
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be used to fit a smooth curve to each sot of contour points. This curve fits 
the interpolated contour points exactly and vai'ies smoothly between them. The 
tension factor may be adjusted to yield suitable curves (a large factor 
corresponds to straight lines between contour points) . The contours which 
results from this technique are, of course, only approximate. Their accuracy 
is limited by the linear interpolating plane procedure and is dependent upon 
the density of the available data. 

PROGRAM DESCRIPTION 
General 

A digital computer program has been developed to implement the contouring 
technique. Program CONTOUR is coded in FORTRAN for a CDC 6600 computer using 
the LRC operating system. Plotting is done using CALCOMP routines and any 
of the graphics post-processors available at LRC. The main program handles 
input and controls primary-level subroutines which contour the data. A brief 
description of the primary-level subroutines is presented in the appendix. 

A listing and flow chart of the main program are given in figures 5 and 6. 

Data transfer is through labeled COMMON and subroutine argument lists. 

Table I describes the COMMON blocks. 

Program Input 

The program is controlled by three things*. 

1. Tape input is used for the data points to be contoured. 

2. Card input defines the program logic. 

3. Field length controls the maximum problem size. 


The data points to be contoured are read from a magnetic tape or disk 
file identified as TAPE 9, The data tape must have been created using the 
LRC output routine RECOUT. The user should prepare TAPE 9 with one record 
for each data point. Each record contains the following: 

X - X-coordinate of system 

Y - Y-coordinate of system 

Z - Value to be contoured 

V - Alternate value to be contoured for same X and Y (if there is 

none, set V = 0) , 

Card input to the program consists of a problem title card and NAMELIST 
input. The problem title card is read with a 3 A10 format and is required 
only once for each run. The NAMELIST input is identified as CASE and each 
successive contour plot requires a separate NAMELIST input. Default values 
and descriptions of each variable of the NAMELIST are contained in Table II. 

If the user desires to use default values, only two variables must be specified 

in the NAMELIST. They are NCL which is the number of contour levels, and CL 

which is a vector containing the actual values of each contour level. NCL is 
presently limited to a maximum of 20. The remaining variables in the NAMELIST 
do not have to be specified unless the user desires to change the default 
values . 

In order to reduce core storage requirements, the program partitions the 
data to be contoured. This partitioning is handled internally and is controlled 
by three variables--IMAX, MAXPTS, and LPCT. IMAX is the maximum number of 
points in each partition. MAXPTS is the total number of points to be plotted 
(equal to the number of records on TAPE 9 if no records are skipped. Note-- 


MAXP'rs s 0 yields a plot of the complete file) . LPCT Is the percentage of 
overlap between partitions in order to yield continuity between partitions. 

Using default values, the program will place a maximum of 1000 data points in 
each partition with a 20 percent overlap, and the complete file will be plotted. 

The field length required for CONTOUR is dependent upon the number of 
data points to be contoured. The program requires 42K octal storage locations. 
To this must be added the larger of 24 times IMAX, or four times MAXPTS plus 
2K octal. Thus, the user may trade-off IMAX and MAXPTS to reduce core 
storage. 

Several options are available to the user. IZPLOT controls the choice 
of which value in the data record is to be contoured (Z or V) . IPLOT 
controls whether contours are drawn by polygonal lines, or by a spline under 
tension (smooth curves). And NDEC and IITNUM control the contour level labels 
whxch may be suppressed if desired. 

To terminate the program, an end of file (EOF) is checked on card 
input. If an' EOF is detected, the plot file is closed and the program is 
terminated. 


Program Output 

Program output consists of user messages and the contour plots. The 
plots are generated by one of the LRC graphics post-processor (e.g., VAR IAN ) . 
Each NAMELIST is printed, along with the total number of points plotted and 
the partition sizes. Data error messages are provided. The error messages 
are contained in the primary-level subroutine descriptions in the appendix. 
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Sainpln .’.•: , *uAv?, 

’ uvir sample plots were gou-valtd o;,d arc- pyo?.f*»tod in ?i p.ures 7 thrash 
10, 'i!iio data used tor tiio plots are typical oi computed rc-sut . r from a 
parametric study of wave conditions. In general, default values of the 
NAMELIST were used and the NAMcLIST required for each plot is given in the 
figures . 

i • 

The figures are identified' as follows: 

figure 7: Ono contour (NCL = 1) of wave height (IZPLOT = 1) smoothed 

with spline under tension (IPLOT = 1) , 

Figure 8: One contour of orbital velocity (IZPLOT = 2 ) smoothed with 

spline under tension. 

Figure 9: Three contours (NCL = 3) of wave height smoothed with spline <* 

under tension. 

Figure 10: One contour of wave heights with polygonal lines (IPLOT = 0). 

For this sample case, there were 3500 data points. The points were 
split into partitions with a maximum size of 1000 points and 20 percent overlap. 

The case required a storage of 121000 octal locations and a run time of 500 
CPU seconds. These values will vary with problem size. 

CONCLUDING REMARKS 

A contouring technique has been presented which can contour randomly 
spaced input data. A digital computer program has been developed at the 
Langley Research Center which implements this technique. The program can 
handle up to 56,000 data points, and provide for up to 20 contour intervals 
for a multiple number of parameters. Contour lines can be either polygonal 
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or smooth, using rpl Jn*' unio 1 - tension. With, relatively minor internal 
modifications, tlr program can handle essentially unl;<';it\,d amounts »£ data 
from a wide rang/; of applicat^or:* . 
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TAOLU I.- DLoCRIiTlON OP ELEMENT iN LABELED « ;)UuW h ? .OCK 


LabUod Common Nine 

tit j Vi “ '•'rjno 

D -..ription 

i \lv\' 

NB 

i.'tpb i* of boundnxy coin s - f the 
tri.uigular grM 


NPTS 

Number of point s of tb i crianguln .* grid 


XSCALO 

X scale factor 


YSCALB 

Y scale factor 


XMIN 

Minimum X value 


YMIN 

Minimum Y value 


XMAX 

Maximum X value 


YMAX 

Maximum Y value 


NL 

Number of lines of the triangular grid 


NT 

Number of triangles of the triangular 
grid 


NI 

Number of interior points of the 
triangular grid 

SIGI 

SIGMA 

Tension factor applied to cubic 
spline 

CONTRL 

CLOVEL 

Contour level value 


NDEC 

Number of digits to the right of 
decimal point on the contour values 


IITN 

Scaled heights of plotted contour 
level digits 

FLAGS 

IPLOT 

Type of smoothing 


ISKIP 

Number of records skipped on the 
input file 
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TABLE I,- CONTINUER 


tghe lad Conaor t Na me yn ti-tblo Name 

I2PLOT 

IMAX 

MAXPTS 

LPCT 

BNDPLT XBMIN 

YBMAX 

XO 
YO 
I PEN 

RAF INDEX(lOl) 

LNG(20) 

BND(2 , 20) 

NBI (2,20) 

IFET(18) 

NBLK 

TMCOM2 LL(4) 


De scripti on 

Plot. typo Flag 

Maximum point" por partition 

Maximum points plotted 

Plot percentage overlap 

.Minimum boundary of the partitioned 
plot 

Maximum boundary of the partitioned 
plot 

Previous plotted X value 

Previous plotted Y value 

Pen position of the plot 

Random Access File index 

Record length for each partition 

Plot boundaries for each partition 

Number of boundary and Interior 
points of each partition 

File Environment Table for the 
Random Access File 

Number of partitions to be plotted 

Communication indices of the line 
array 



TABLE II,- DESCRIPTION OF ELEMENTS IN NAMELIST CASE 


Variaolc Name 

Default Value 

Description 

NCL 

- 

Number of rontour levels (maximum 
of 20) 

CL 

- 

Actual values of contour levels 
(dimensioned 20) 

I PLOT 

1 

= 0 yields polygonal lino 
contour 

= 1 yields smooth contours 
(splino with tension) 

IS, UP 

0 

Number of records to skip 
between input data points 

IZPLOT 

1 

= 1 plot of z-values 
= 2 plot of v-values 

IMAX 

1000 

Maximum number of points to be 
plotted per partition 

MAXPTS 

0 

Maximum total points to be 
plotted 

= 0 Plot the complete file 

NDEC 

- 1 

Number of digits to the right 
of decimal point on the contour 
= - 1 Suppresses decimal 

HTNUM 

.0,15 

Height of plotted contour level 
digits (inches) 

= 0 Suppresses contour numbers 

LPCT 

20 

Percent overlap of the plot 
partitions 
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TABLE III - DESCRIPTION OF SECONDARY SUBROUTINES 


Subroutine Name Description 


REQFL 

Finds out present program field 
length or requests a larger or 
smaller field length. 

INITFET 

Initializes a File Environment Table. 

ADDFILE 

Adds a file name to the communications 
area.. 

SUBFILE 

Subtracts a file name from the 
communications area. 

LTEST 

Computes the quadrant of a point. 

STEST 

Tests two triangles to maximize the 
minimum angles. 

KURV1 

Computes parameters necessary for a 
spline under tension for an open 
curve . 

KURV2 

Computes the points for a spline 
under tension for an open curve. 

i 

KURVP1 

Computes parameters necessary for a 
spline under tension for a closed 
curve , 

KURVP2 

Computes the points for a spline 
under tension for a closed curve. 

FRSTPT 

Moves plotter to the first point 
of a line. 

VECTOR 

Moves plotter to successive points 
of a line. 

POINT 

Plots a single point. 
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TABLE IV - DESCRIPTION OF SYSTEM ROUTINES 


Routine Name 
RECIN 

OPBNMS 

WRITMS 

READMS 

EVICT 

LOCF 


Description _____ 

Reads binary records written by the 
LRC routine RECOUT. 

Initializes a random access file. 

Writes a random access file. 

Reads a random access file. 

Deletes a file from the system and 
releases the disk space. 

Defines the address of a word in core. 
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PROGRAM CONTOUR! |NPUT*1 .OUTPUT*! . TAPES* INPUT « TAPE6»OUTPUT » 
DIMENSION C00El3l.Ct.IZ0) 

COMMON /XTXRX I / NO . NPTS • XSCAlE . VSCALE.XM I N . VM IN. XMAX . YMAX . NL .NT <NI 

COMMON /SIGI/ SlT.MA 

COMMON /CONTRL/ CLEVEL .NDEC . HTN 

COMMON ✓FLAGS,' IPLOT. ISKIP. IZPLOT. IMAX. MAXPT5.LPCT 
COMMON /HNDPLT/ XPMlN.XOMAX.XO. VO. |‘PEN 

COMMON ✓RAF/ INDEXtlOt ) . LNGIZO I .QN0I2.30 ) .NO! (Z.ZOt.lFET! IBI.N0LK 
COMMON A(J) 

DYNAMIC STORAGE IS USED FOR STORAGE OF XYZ POINTER. LINE AND 
TRIANGLE ARRAYS' SO PROGRAM FIELD LENGTH DEFINES CORE STORAGE 


FL ■ A2* OCTAL * MAX I MUM I Z4» IMAX . A»MAXPTS ) DECIMAL 
NAMELIST INPUT 

NCL -NUMBER OF CONTOUR LEVELS (MAX ■ SO I 

CL -ACTUAL VALUES OF CONTOUR LEVELS 

I PLOT -*0 YIELDS POLYGONAL LINE CONTOUR 

*1 YIELDS SMOOTH CONTOURS (SPLINE WITH TENSION! 
ISKIP -NUMBER OF RECORDS TO SKIP BETWEEN ACTUAL DATA PT5 
I2PL0T -»l WAVE HEIGHT PLOT 
*2 ORBITAL VEL PLOT 

IMAX -MAXIMUM NUMBER OF POINTS TO BE PLOTTED IN ONE BLK 
MAXPTS -MAXIMUM TOTAL POINTS TO BE PLOTTED 
*0 PLOT COMPLETE FILE 

NOEC -NUMBER OF DIGITS TO RIGHT OF DECIMAL POINT ON 
CONTOUR VALUES (-1 SUPPRESSES DECIMAL) 

HTNUM -HEIGHT OF PLOTTED CONTOUR LEVEL DIGITS (INCHES) 

=0 SUPPRESSES CONTOUR NUMBERS 
LPCT -PERCENT OVERLAP OF THE PARTITIONS 


* 


« 


* 


< 


RANDOM ACCESS FILE STRUCTLWE 
RECORD I FIRST BLOCK 


RECORD 

a 

FIRST BLOCK V 

RECORD 

3 

FIRST BLOCK Z - WAVE HEIGHT 

RECORD 

4 

FIRST BLOCK V - ORBITAL VEL 

RECORD 

5 

SECONO BLOCK X 

RECORD 

ft 

• 

SECONO BLOCK Y 

RECORD 

• 

• 

4»NBLK-3 

LAST BLOCK X 

RECORD 

4*N0LK-2 

LAST BLOCK Y 

RECORD 

4 #N0LK* 1 

LAST BLOCK z - WAVE HEIGHT 

RECORb 

4 •NBLK 

LAST BLOCK V - ORBITAL VEL 

RECORD 

4*NBLK4t 

FIRST BLOCK LINE AND 1 TR 1 ARRAYS 

RECORO 

4«NBLK*2 

• 

SECOND BLOCK LINE AND ITRJ ARRAYS 

RECORD 

• 

• 

5#NBLK 

.ASi BLOCK LINE AND ITRI ARRAYS 


PROGRAM is STRUCTURED FOR A MAXIMUM OF TWENTY PARTITIONS 
WHERE NBLK * I100»(MAXPTS— IMAXI !✓( I I OO-LPC T I • I MAX )♦ 1 

NAMELIST /CASE/ NCL. CL. IPLOT. ISKIP, IZPLOT. IMAX. MAXPTS. NDEC. HTNUM 
I .LPCT 

INITIALIZE NAMELIST 
.IMAX * 1000 
MAXPTS * 0 
NDEC * -I 
ISKIP * 0 
IZPLOT • l 
IPLOT * 1 
HTNUM * 0,15 
LPCT * SO 


INITIALIZE PROGRAM 
SIGMA ■ -3 

sc * ao.o 

XD » SC»S1 ,0 
YD * SC#I0.0 
TB * SC 

ts « a, o/sc 

HT * SC*0,3 
XS = SC*3.0 

ys - sc*a.o 

DV « 10 O/SC 

(USED * 1000000 

01 V " 1 .0E-10 

R » UR AN ( ,31 A 1 $9 1 9S037 A ) 

1 FLAG = O 




i 

j 


INITIALIZE THE PLOTTER 




Figure 5. 


Listing of main 


program 
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CALL PSEUDO 
CALL LEROY 


PROBLEM TITLE 
READ IS. EDO) CODE 
FORMAT I DA | 0 ) 
CONTINUE • 


READ NAMELIST 
READ IS i CASE ) 

IF (EOF, SI 20,30 
CONTINUE 


CLOSE PLOT 
CALL CALPLT 10,0*099 I 
STOP 

CONTINUE 

WR I TE ( 6 , CASE I 

HTN * SC«HTNUM 


SET UP AXES 

CALL -ALPLTIUO.I ,0.-3) 

CALL NOTATEI-AO,0,0,HT, CODE, 90. 0,30 I 
CALL AXESIO.O.O.XD.O.OV.TR.TS.IM .HT.-ll 
IF [ IZPL0T.EQ.2 I GO TO An 

CALL AXE5(0,0,gO,0,YD.O,DU.TB,TS,20HWAVE HEIGHT (METERS ) , HT.20 I 

GO TO SO 

CONTINUE 

CALL AXESI0.0,90.0.YD<0.0V.TB,1 S.20H0RBITAL VEL ( CM/SEC I »HT ,20) 
CONTINUE 

CALL GRIDI0,0,XS,YS,7,3) 

IF I IFLAG.GT.O I GO TO 70 
1FLAG * 1 


READ THE INITIAL POINT AND WRITE ON RAF 
CALL READ IN 


LOOP ON THE NUMBER OF BLOCKS 
00 60 1*1 , NBLK 
IREC * 4« t t-t 1 + 1 
N * LNG ( I ) 


LOCATE ALL ARRAYS IN CORE 
LOCX ■ t 
LOCY ■ LOCX+N 
LOCP ■> LCCY+N 
LOCL ■ LOCPfN+2 


READ -x and Y from RAF 
CALL READMS ( 9 • A (LOCX ) ,N, I REC I 
IREC * IREC+I 

CALL READMS ( 9 . A (LOCY ) . N, I REC I 


COMPUTE THE TRIANGULAR GRID 
CALL TRI ANG( A (LOCX ) .A (LOCY ) ,N, A (LOCP ) , 1SUP, I USED 1 


WRITE THE LINE AND TRIANGLE DATA TO RAF 
I RFC = «*NBLK+t 
NR I < 1 , I I » N9 
NB t I 2 , I ) * N I 
NLC = 2»ND+3»(N-NB-I > 

NTC “ NB+2»(N-NB-1 | 

LENGTH '» 5»(NLC+i |+3*(NTCfl > 

CALL WR I TMS I 9 . A ( LOCL ) • LENGTH ■ I REC 1 

CONT I NUE 

CONTINUE 


LOOP ON THE BLOCKS 
DO 100 !°t,NBLK 
IREC ■ «• ( l-l *+l 


THE LENGTH FOR THE XYZ RECORDS 
N = LNG ( I > 


LOCATE THE ARRAYS IN CORE 
LOCX » I 
LOCY « LOCX+N 
LOCZ ■ LOCY+N 
LOCL ■ LOCZ+N 


READ THF X ARRAY 
CALL READMR( 9, A I LOCX ) ,N, 1 REC ) 
I RFC » IRFC+1 


READ THE V array 
CALL HEADMS(9,A(L0CYI*N.IREC) 
IREC - IREC+IZPLOT 


Figure 5. Continued. 


2BSS 


* READ THF l 0« V ARRAY 
CAUL READM3I9.A ILOCZ ) >N< |REC| 

• 

• PLOT! INC BOUNDARIES FOR THE PARTITION 
XOMIN « ONO I I . I | 

XBMAX • BNDI2. I I 

* 

• DEFINE THE NUMBER OF BOUNDARY POINTS NQ 

• INTERIOR POINTS N| 

• GRID LINES Nl 

• GRID TRIANGLES NT 

NR » NQ I 1 1 • I ) 

Nl « N0I (2. I ) 

NL ■ 2«ND+3»(NI-1 | 

NT « NB+2MNI-I I 
NLC ■ 2»NB+3« IN-NB-I | 

NTC ■ NQ*2» (N-NB-I ) 

LENGTH * S»INLC+I 1+3»INTC+I) 

• 

* READ the LINE and triangle, arrays 

|REC « A*N9lK+| 

call READMS 1 9 . A I LOCL ). LENGTH . I REC I 
IF ( I PLOT .GT ■ 0 ) GO TO VO 

• 

* POLYGONAL LINE CONTOUR PLOTS 
2 SPEC » I.OE+9 

call TRIOROIAfLOCZI .N. AiLOCL) .ZSPEC) 

DO BO J»| >NCL 
CLEVEL * CLIJ) 

CALL SCAN ( A I LOCK ) I A I LOCY I • A (LOCZ I .N. A ILOCL ) . CL ( J I I 
BO CONTINUE 
GO TO 100 
90 CONTINUE 

# 

» SMOOTH CONTOUR PLOTS (SPLINE WITH TENSION) 

I .OCT * LOCL+5* (NLC+1 ) 

CALL CONTRI < A ILOCX ) . A (LOCY t . A ILOCZ ) . A (LOCL > . A (LOCT ) .CL .NCL.MODE I 
[F (MODE.EO.l) GO TO 100 
WRITE (6. 20 1 | MODE 

201 FORMAT I# — -CONTRI ERROR. MODE ■» 1 2 ) 

100 CONTINUE 

• 

» NEXT FRAME 

CALL NFRAME (20.0.0 ) 

GO TO 10 

end 




Pa GB In 

QUAifg 


Figure 5, Concluded 
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CLOSE PLOT 
CALL CALPL^ 


CONSTRUCT 
CASE AXIS 


NEXT FRAME 
CALL NFRAME 


FIRST 

PASS 


READ DATA SET 
AND CREATE RAF 
CALL READIN 
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Figure 6. Flow chart of main program 

















Namelist input: 


$CASF. 

NCL*= 1 , 

CL-O. 1E+01, 
IPLOT*= 1 , 

ISKIP=0 , 

I ZPLOT= 1 , 

IMAX* 1 00 0 # 
MAXPTS*0 , 
NDEO-1 , 
HTNUM=0. 15E+00, 
SEND 


Figure 7 . One contour of wave height smoothed 
with spline under tension. 19 





)ne contour of orbital velocity smoothed 
with spline under tension. 


Namelist input: 


Figure 8. 

2 " 


$CASE 
NCL=1, 

CL-0.15E+02, 

I PLOT-1, 

1SK IP=0 , 

. ZPL0T=2 , 

IMAX* 1000 , 
MAXPTS=0 , 
NDEC— 1 , 
HTNUM*0 . 15E+00 
SEND 







^ • 

V 





Numeliat L.put: 


$CASE 
NCL=»3 , 

CL-O. lE*01 f 0.2E+01 ,0. 3E+01, 
I PLOT- 1, 

ISKIP-0 , 

IZPLOT-1, 

IMAX-1000 , 

MAXPTS=0, 

ndeo-1, 

HTNUM=0 . 15E+00 , 

$END 


I 


Figure 9. Three contours of wave height smoothed 
with spline under tension. 
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Namelist Input: 


$CASE 

NCL-1, 

CL-0 . 1E+01, 

I PLOT- 0 , 
ISKIP-0 , 
IZPLOT- 1 , 
IMAX-1000, 
MAXPTS-0 , 
NDEC-- 1 , 
HTNUM-0 . 15E+00 
$END 


/>» 


'.r n 
+QC v 


Or 


An 


1 

Figure 10. One contour of wave heights with 
polygonal lines. 
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APPENDIX 

Description of Primary - Levol Subroutines 

The primary-level subroutines are those called by the 
main program and major modules of the respective basic steps 
of the program. Brief descriptions of the secondary subroutines 
and system routines are contained in Tables III and xV. 

The following is a description and use of each of the 
primary subroutines. 


Noma: 


RGADIN 


Description: 

Use : 

Labeled Common: 
Procedure : 


External Routines’: 

Name: 

Description: 


Subroutine READIN is a driver routine 
for the contour data input and storage 
routines. Dynamic storage allocation is 
used and the program field length is 
defined accordingly. 

CALL READIN 

/RAF/JNDEX(101), LNG(20), BND(2,20) , 

NBI (2,20) , IPET(18), NBLK 

/ FLAGS/ IPLOT , ISKI P , I ZPLOT , IMAX , MAXPTS , LPCT 

An input file is created using dynamic 
storage for the FET and buffer. The 
data points are then road into the 
available core (routine GETXYZ) . The 
minimum and maximum values and the scale 
factors of the X and Y coordinate arrays 
are computed (routine MINMAX) . The arrays 
are ordered with a major sort on increasing 
X and a minor sort on decreasing Y, with 
duplicate values deleted from the set (rou- 
tine SHELL). The input file is released 
and a random access file, RAF, is created. 

The input data are stored in partitioned 
form on RAF (routine PUTXYZ). The 
field length is then reduced to the required 
level to execute the program. 

GETXYZ, MINMAX, SHELL, PUTXYZ, REQFL, ADDFILE, 
INITFET, SUBFILE, EVICT, OPENMS, LOCF 


GETXYZ 


Subroutine GETXYZ reads the contour input 
data from a disk or a magnetic tape. The 
file read must be generated by the LRC 
output routine RECOUT. 
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CALL GETXYZ (X, Y, Z, V) 


Use : 


Labeled Common: 
Procedure : 


External Routines : 

Name : 

Description: 

Use: 

Labeled Common: 


X - An array of X-coordinatos of the 
contour plot. 

Y - An array of Y-coordinatos of the 

contour plot. 

Z - An array of wave heights. 

V - An array of orbital velocities, 

/FLAGS/ I PLOT, ISKIP , IZPLQT , IMAX ,MAXPTS , LPCT 

Entries are made into the X,Y,Z and V 
arrays until the limit of the points 
to be read is exceeded or an end of file 
is detected. The capability of skipping 
a desired number of records is provided. 

If on end of file is detected before the 
limit of points is exceeded, the limit is 
redefined to the number of points on the 
file. 

RECIN 


MINMAX 


Subroutine MINMAX identifies the minimum 
and maximum values of two arrays, and 
computes a scale factor for each array. 

CALL MINMAX (X,Y) 

X - Array of values 

Y - Array of values 

/XTXRXI/NB ,NPTS,ISCALE,YSCALE,XMIN,YMIN 
XMAX,YMAX,NL,NT,NI 


/FLAGS/IPLOT, ISKIP, IZPLOT, IMAX, MAXPTS, LPCT 



Procedure : 


Error Conditions: 


Name : 


Description : 


Use : 


Procedure : 


The two arrays, X and Y, are searched 
for the minimum and maximum values of each. 
The differences between the respective 
minimum and maximum values are used to 
compute the scale factors for the X and 
Y arrays. 

If the difference between a minimum and a 
maximum value is zero, a message is given 
and the program is terminated. 


SMELL 


Subroutine SHELL uses the shell sorting 
technique to perform a major and minor 
sort, disregarding duplicate points. 

CALL SHELL (X,Y,Z,V,N) 

X - Array of the major sort. 

Y - Array of the minor sort. 

Z - Array sorted with respect to X and Y. 

V - Array sorted with respect to X and Y. 

N - Number of elements per array sorted. 

A shell sorting technique is used to sort 
the four arrays with a major sort on 
increasing X and a minor sort on decreas- 
ing Y. If duplicate X-Y points are 
detected one set of points is deleted 
and the sort is reinitialized. 
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Name : 


PUTXYZ 


Description: 


Use: 


Labeled Common: 


Procedure : 


Error Condition: 


External Routines : 


Subroutine PUTXYZ defines partitions for 
plotting and stores each partition on a 
random access file. 

CALL PUTXYZ (X,Y,Z,V) 

X - Array to be partitioned. 

Y - Array to be partitioned. 

Z - Array to be partitioned. 

V - Array to be partitioned. 

/ F L AGS / 1 P LOT, , I S KI P , I ZP LOT , IM AX , MAXPTS , LP CT 

/RAF/ INDEX (101) ,LNG(20) ,BND(2,20) , 

NBI (2,20) , IFET(18) ,NBLK 

The number of partitions are defined 
such that there will be an overlap of 
points between partitions. The partition 
sizes are computed to have an equal number 
of points per partition with an input over- 
lap. For each partition four records are 
written for the respective arrays . The record 
lengths and plot boundaries of the overlaps 
are stored in core. 

If the number of partitions exceed the 
program limit of twenty, the program is 
terminated. 

WRITMS 


Name : 


TRIANG 


Description : 


Subroutine TRIANG is the driver routine for 
the generation of the triangular grid. 


Use: 


CALL TRIANG (X,Y,N,1TRI>ISUP,1USED) 


Labeled Common: 
Procedure : 


External Routines : 


Name : 


Description : 


Use: 


X - Array of X-coordinated points for the 
triangular grid. 

Y - Array of Y-coordinate points foT the 
triangular grid. 

N - Number of points of the X and Y arrays. 

ITRI - Work array containing the triangulariza- 
.information. 

ISUP - Dimension of the ITRI array 

IUSED - Actual number of data locations 
required for the ITRI array. 

/XTXRXI/NB ,NPTS , XSCALE , YSCALE , XMIN 
YMIN,XMAX,YMAX,NL,NT,NI 

The boundary points of the triangular 
grid are defined (routine CONHUL) . An 
initial grid is generated (routine TMESH2) 
and then an iteration on the grid is 
performed for improvement (routine TMESH3) . 

CONHUL, TMESH2,TMESH3 


CONHUL 

Subroutine CONHUL identifies the boundary 
points for the triangular grid such that 
the boundary points form a convex hull. 

CALL CONHUL (X,Y,NPTS,IP,IP2 ,M0DE) 

X - Array of X-coordinate points. 

Y - Array of Y-coordinate points. 

NPTS - Number of points in the X-Y arrays . 
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Labeled Common: 
Procedure : 

Error Conditions : 


Name : 

Description : 

Use : 


IP - A pointer array of length NPTS. 

IP2 - A dummy array. 

MODE - An error flag equal to one if the 
routine terminates normally. 

/XTXRXI/I 1 ,NP , XSCALE , YSCALE , XMIN 
YMIN,XMAX,YMAX,NL,NT,NI 

The IP array is initialized pointing to the 
sorted X-Y arrays. The minimum point is 
considered the first point of the boundary. 
Then the successive points are checked by 
using the cross product, such that the 
boundary is defined in a counter-clockwise 
order. The pointer array is swapped so that 
the first entries contain the indices of the 
boundary points . 

MODE = 2 Less than four points . The 

algorithm requires at least four 
points , three boundary and one 
interior point . 

MODE - 5 No interior points . 

MODE - 6 All values of the X array equal 
the maximum value. 


TMESH2 


Subroutine TMES1I2 computes an initial 
triangular grid for the contour plotting. 

CALL TMESI12 (X, Y, NPTS, IP, LINE, ITRI, MODE) 

X - Array of X-coordinate points for the 
triangular grid. 

Y - Array of Y-coordinate points for the 
triangular grid. 


Labeled Common: 
Procedure : 

Error Conditions: 

External Routines 

Name : 

Description : 

Use: 


NPTS - Number of points in the X and Y arrays. 

IP - Pointer array for the X and Y arrays. 

LINE - Array of descriptors of each line 

ITRI - Array of descriptors of each triangle 
of tiie triangular grid with three 
entries per triangle. 

MODE - An error flag equal to one if the 
routine* terminates normally. 

/XTXRXI/I 1 , NP , XSCALE , YSCALE , XMIN , YMIN , 
XMAX.YMAX, NL, NT, NI 

A centroid of the interior points is computed 
and the nearest point is chosen to initialize 
the triangular grid. Lines and triangles are 
defined from the boundary and the centroid 
points. The remaining interior points are 
introduced creating a complete triangular grid. 

MODE » 6 A point is not contained in any 

triangle and is omitted from the set . 

MODE =7 Bad data is encountered and the 
routine is exited. 

MODE =8 A point not included in the boundary 
set but now lie*; on the boundary. 

The point is omitted from the set. 

STEST, LTEST 


TMESH3 


Subroutine TMESH3 performs an iterative 
improvement of the triangulation . 

CALL TMESH3 (X,Y, LINE, ITRI) 

X - Array of X-coordinate points for the 
triangular grid. 
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Labeled Common: 


Procedure : 


External Routine: 


Name : 


Description: 


Use: 


Labeled Common: 


Y - Array of Y-coordinate points for the 
triangular grid, 

LINE - Array of descriptors of each triangle 
of the triangular grid with three 
entries per triangle. 

/XTXRXI/NB ,NPTS , XSCALE , YSCALE ,XMIN , YMIN 
XMAX,YMAX,NL,NT,NI 

/TMC0M2/LL(4) 

An iteration is performed to improve 
the triangulation by interchanging 
lines of the triangle so that the 
smallest angles of the triangles: are 
maximized. A limit of three times the 
number of lines is imposed on the 
iteration count. 

STEST 


TRIORD 


Subroutine TRIORD reorders the line indices 
with respect to the Z vector, such that 
moving from the initial to the terminal 
point on a line involves an increase 
in the associated Z values. 

CALL Til ORD (Z,N,ITRI,ZSPEC) 

Z - Array associated with the line and 
triangle indices. 

N - Number of values in the Z array. 

ITRI - Array containing the line and 
triangle information. 

ZSPEC - A value of the Z array which is to be 
ignored. 

/XTXRXI/NB , NPTS , XSCALE , YSCALE , XMIN , YMIN , 
XMAX,YMAX,NL,NT,NI 


Procedure : 


Name : 


Description : 


Uso: 


Labeled Common: 
Procedure : 


External Routines : 


For each of the linec » t:he associated Z 
values of the end points are reordered 
such that the terminal Z value is greater 
than the initial Z value. Sign changes 
are made on the corresponding triangle 
data if there is a reordering. 


SCAN 


Subroutine SCAN produces the orderod 
sequence of X-Y pairs of points on a 
given contour. 

CALL SCAN (X,Y,Z,N,ITRI,CL) 

X - Array of X coordinates of the surface 
data. 

Y - Array of Y coordinates of the surface 
data. 

Z - Array of Z coordinates of the surface 
data. 

N - Number of data points . 

ITRI - Array containing the line and triangle 
information. 

CL - Value of the contour. 

/XTXRXI/NB ,NPTS , XSCALE , YSCALE ,XMIN 
YMIN,XMAX,YMAX,NL,NT,NI 

For each line, a check is made for an 
intersection. The contour is continued 
by checking the lines of adjacent triangles 
until the contour is enclosed or intersects 
a boundary line. If the boundary is 
intersected, the contour is continued 
backwards from the first point until the 
second boundary is encountered. For 
each line crossed, a linear interpolation 
between the end points is made and commands 
are given to plot the computed point. 

FRSTPT, VECTOR 
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Namo: 


CONTRI 


l 

' Description: 

Use : 


Labeled Commpn: 
Internal Data: 

Procedure: 


Subroutine CONTRI checks the triangular 
grid data and initializes the contour 
arrays for the spline with tension curve 
fitting. 

CALL CONTRI (X,Y,Z,LINE,ITRI,CTRVAL,NC,MODE) 

X - Array of X-coordinate points . 

Y - Array of Y-coordinate points. 

Z - Array of Z-coordinate points, 

LINE - Array of descriptors of each line of 
the triangular grid with five entries 
per line. 

ITRI - Array of descriptors of each triangle 
of the triangular grid with three 
entries per line. 

CTRVAL - Array of contour levels to be plotted. 

NC - Number of contour levels. 

MODE - An error flag equal to one if the 
routine terminates normally. 

/XTXRXI/I 1 , NP , XSCALE , YSCALE , XMIN , YMIN , 
XMAX,YMAX,NL,NT,NI 

/CONTROL/CLEVEL , NDEC ,HTN 

XPS - Array of interpolated X values 
associated with a contour level, 
dimensioned 500. 

YPS - Array of interpolated Y values 
associated with a contour level, 
dimensioned 500. 

The line descriptor array is reordered 
such that the end points of the line are 
increasing with respect to the Z values. 

For each contour level, the line descriptor 
array is searched for an intersection, 
defining a starting point of a contour, 

A neighboring triangle is used to continue 
the contour until a boundary is reached 
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Error Conditions : 


External Routine: 


Name: 


Description : 


Use : 


Internal Data: 


or the contour is closed. When a line is 
used for a contour, it is flagged so it 
will not be used a second time. Arrays 
are generated from the interpolated X-Y 
values of the end points of the lines 
crossed to be used in the splint' with 
tension curve fitting routines (routine 
CURVL) . 

MODE * 7 Bad data structure. A line 
does not have an associated 
triangle. 

MODE = 8 Bad data structure. A contour 
level crossing one side of a 
triangle does not cross another. 

MODE = 9 Too many lines crossed by a 

contour. The XPS and YPS arrays 
are insufficiently dimensioned. 

CURVL 


CURVL 


Subroutine CURVL determines if the curve 
is closed or open and calls the appropriate 
spline routines. 

CALL CURVL (X,Y,NPPTS,NNPTS,FAC) 

X - Array of X-coordinate points to be 
curve fitted. 

Y - Array of Y-coordinate points to be 
curve fitted, 

NPPTS - Number of points of the contour 

before a boundary line is encountered 
in the arrays. 

NNPTS - Number of points of the contour 

after a boundary line is encountered 
in the arrays . 

FAC - The magnitude of the plot. 

XP - Array for the storage of curvature 
information for the given nodes, 
dimensioned SOD. 
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SCRAT - Temporary storage array, dimensioned 1000. 


labeled Common: 
Procedure 


External Routines; 


/SIG1/SIGMA 

If the arrays define a closed contour, 
the closed curvo routines are executed 
(KURVP1 and KURVP2) . If the contour is 
open, the arrays are reordered such that 
the first and last points are the end 
points , and the open and curve routines 
are executed CKURV1 and KURV2) . When 
a point is defined the plotting routines 
are executed, 

KURV1, KURV2 , KURVP 1 , KURVP2, POINT, I’RSTPT, 
VECTOR. 
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